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Abstract 


Application  of  the  Finite  Element  Method 
to  Random  Rough  Surface  Scattering 
with  Neumann  Boundary  Conditions 

by  Kevin  Krause 


Chairperson  of  Supervisory  Committee: 


Abstract: 


Professor  Leung  Tsang 
Department  of  Electrical  Engineering 


Scattering  from  a  one;dimensional  rough  surface  with  Gaussian  roughness  spec¬ 
trum  is  analyzed  using  js.  finite  element  formulation.  The  method  is  applied  to  Monte 
Carlo  simulations  satisfying  Neumann  boundary  conditions.  Finite  element  results 
are  compared  with  Results  obtained  by  solving  an  integral  equation.  Convergence  of 
the  method  is  verified  by  varying  the  number  of  nodal  points  in  the  first  order,  trian¬ 
gular  mesh.  Results  are  in  excellent  agreement  with  tapered  wave  integral  equation 
solutions  for  large  surface  length  after  averaging  over  realizations.  Finite  element  ad¬ 
vantages  in  •CPU' time  and  memory  storage  are  presented  for  the  examples  discussed. 
Comparisons  with  the  Kirchoff  approximation  and  small  perturbation  theory  within 
their  respective  regions  of  validity  are  also  presented.  Additionally,  analysis  of  the  ef¬ 
fects  of  decreasing  surface  length  jn  incoherent  scattering  results  of  the  finite  element 
method  is  accomplished  to  investigate  the  method’s  likely  advantages  for  large-scale 
scattering  problems.  Numerical  Results  of  scattering  are  represented  in  terms  of  the 
normaiized  bistatic  scattering  cross  section. 
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Chapter  1 

INTRODUCTION 


Scattering  of  electromagnetic  waves  from  rough  surfaces  has  been  extensively  stud¬ 
ied  using  the  Rayleigh-Rice  perturbation  theory  and  the  Kirchoff  approximation 
[2,6,16,17].  These  classical  approaches  assume  plane  wave  incidence  upon  arbitrar¬ 
ily  long  surfaces  and  utilize  statistical  properties  in  calculating  ensemble  averages. 
Each  is,  however,  limited  in  its  domain  of  validity  which  depends  upon  the  surface 
characteristics. 

With  the  advent  of  modern  computers,  increased  interest  in  Monte  Carlo  simula¬ 
tions  of  scattering  by  random  rough  surfaces  has  evolved.  In  Monte  Carlo  simulations, 
ensemble  averages  are  calculated  by  averaging  scattered  field  intensity  over  hundreds 
of  surface  realizations.  Computers  are  well  suited  for  the  task  of  averaging  repetitive 
calculations,  nonetheless,  it  is  important  to  minimize  the  computation  time  for  each 
realization. 

The  most  common  Monte  Carlo  method  for  an  “exact”  numerical  result  has  been 
the  solution  of  the  tapered  wave  integral  equation  by  the  method  of  moments  [1,3, 
5,16,17],  A  plane  wave  is  tapered  to  avoid  edge  effects  from  a  finite  surface  using  a 
Gaussian  taper  function  which  results  in  an  incident  wave  consisting  of  an  angular 
spectrum  of  plane  waves  about  an  average  incident  angle.  Since  the  width  of  the 
angular  spectrum  is  inversely  proportional  to  surface  length,  large  surface  length 
is  required  to  better  approximate  plane  wave  incidence.  The  disadvantage  of  the 
integral  equation  approach  is  that  it  requires  solution  of  a  full,  complex  matrix, 
and  for  large  scale  rough  surface  problems,  applicability  of  the  method  is  limited  by 
available  computer  memory  and  CPU  time. 

In  order  to  address  these  problems,  application  of  the  finite  element  method  to 
Monte  Carlo  simulations  of  rough  surface  scattering  has  recently  been  investigated  for 
one-dimensional  Gaussian  rough  surfaces  with  periodic,  Dirichlet  boundary  conditions 
(8,9j.  In  this  thesis,  existing  finite  element  code  for  Dirichlet  boundary  conditions  is 
converted  and  applied  to  Monte  Carlo  simulations  of  one-dimensional  Gaussian  rough 
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surfaces  with  periodic,  Neumann  boundary  conditions.  Accuracy  of  toe  results  is  veri¬ 
fied  by  comparison  with  “exact”  results  determined  through  the  use  of  a  tapered  wave 
integral  equation  (TWIE)  solution  with  large  surface  length.  Two  propagating  Flo- 
quet  modes  exist  for  each  wavelength  of  surface  length  in  the  finite  element  approach. 
Since  evanescent  waves  rapidly  decay  above  the  rough  surface,  a  total  number  of  only 
four  to  six  evanescent  modes  are  needed.  In  contrast,  discretization  in  the  integral 
equation  approach  requires  ten  points  per  wavelength.  Hence,  for  a  surface  with 
Lj A  =  30,  the  TWIE  approach  requires  300  basis  functions  while  66  modes  are  suf¬ 
ficient  in  the  finite  element  approach  [9j.  Additionally,  a  full,  complex  matrix  must 
be  solved  for  the  TWIE  approach  while  a  sparse,  real  matrix  must  be  solved  for  the 
FEM  approach.  Thus,  TWIE  generally  requires  significantly  more  computer  CPU 
time  and  memory  storage.  Surfaces  are  generated  using  Monte  Carlo  surface  genera¬ 
tion  routines  with  Gaussian  roughness  spectrum  and  specified  statistics.  Additional 
comparisons  with  both  the  Rayleigh-Rice  theory  and  the  Kirchoff  approximation  in 
their  respective  domains  of  validity  are  also  made.  Finally,  tests  are  conducted  to 
analyze  the  effects  of  decreasing  surface  length  on  incoherent  scattering  results  of  the 
finite  element  method.  This  is  an  important  investigation  since  incoherent  scattering 
is  the  measured  quantity  in  monostatic  and  bistatic  scattering,  and  analysis  of  smaller 
surface  lengths  requires  much  less  CPU  time  and  memory  storage.  The  scattered  field 
quantity  calculated  is  the  normalized  bistatic  scattering  cross  section. 


Chapter  2 

FINITE  ELEMENT  METHOD 


In  the  finite  element  approach, .scattering  from  an  incident  plane  wave  impinging 
upon  a  one-dimensional,  random  rough  surface  with  Gaussian  statistics  is  considered. 
Periodic  boundary  conditions  are  utilized  to  truncate  the  finite  element  mesh  and  to 
discretize  the  plane  wave  spectrum  [2].  In  the  formulation,  the  scattered  field  in  the 
region  above  the  maximum  height  of  the  rough  surface  is  expressed  in  terms  of  Floquet 
modes.  The  finite  element  method  is  utilized  to  solve  the  wave  equation  below  the 
maximum  height.  The  attraction  of  the  finite  element  method  is  the  banded  nature  of 
the  resulting  matrix  equation  and  reduced  memory  storage  requirements  over  those  of 
the  integral  equation  approach.  Formulation  of  the  finite  element  method  is  described 
in  detail  in  reference  (9]  and  will  only  be  briefly  summarized  here. 

2.1  Formulation 

The  geometry  of  scattering  from  a  one-dimensional,  periodic  rough  surface  (two- 
dimensional  scattering)  in  the  xz  plane  is  considered  as  illustrated  in  Figure  2.1.  The 
surface  has  a  Gaussian  roughness  spectrum,  \V{K)  =  (hil/2\/z)exp(—K3l2/4),  and 
the  random  surface  height  is  described  by  z  =  f(x)  with  rms  surface  height  h  and 
correlation  length  /.  The  sample  surface  length  L  is  extended  periodically  with  period 
h  in  order  that  periodic  boundary  conditions  may  be  used.  A  plane  wave  is  incident 
on  the  surface  with  angle  0,  from  the  vertical,  and  Neumann  boundary  conditions 
are  imposed  on  the  surface.  Above  the  surface,  the  xz  plane  is  divided  into  regions  I 
and  I!,  the  homogeneous  and  inhomogeneous  regions  respectively.  Region  I  consists 
of  z  >  d  and  region  II  consists  of  d  >  z  >  /( i),  with  d  chosen  to  be  larger  that  the 
maximum  height  of  the  rough  surface. 

Ail  field  quantities  are  assumed  to  have  e~iu‘  time  dependence,  and  the  incident 
field  is  a  plane  wave  of  the  following  form: 

ik,tz-ik{gz 


tpi(x,z)  =  e' 


(2.1) 
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where  k,x  =  ks'mO,,  k„  --  k  cos  0,  and  k  =  tv^/ioCo  =  2~/A.  The  incident  field  here 
differs  from  that  used  in  the  formulation  of  reference  [9]  by  the  factor  e~,k"d.  This 
difference  insures  the  same  incident  field  for  each  realization  despite  variation  in- the 
parameter  d — a  necessary  condition  for  the  calculation  of  the  incoherent  scattering 
cross  section.  Use  of  periodic  boundary  conditions  results  in  a  discretization  of  the 
scattered  plane  wave  spectrum  into  Floquet  modes  [7].  The  scattered  fields  in  region 
I  consist  of  upward  going  waves  only  and  can  be  expressed  as: 

<!>l(x,z)  =  f]  (2.2) 

771  = —CO 


with 


kxm  — 


kix  + 


2:rm 

~T 


V-k\m>  0 
&.-**>  0 


m  =  0,±1,±2 .  Modal  field  amplitudes  Bm  are  unknowns  which  are  to  be 

determined.  From  Floquet ’s  theorem,  +  L,z)  =  e'k[‘Ltp,(x,z),  and  with  dis¬ 
crete  incident  angle  index  m0,  k,x  can  be  discretized  by  k,x  =  2xm0/L  so  that 
+  L,z)  =  Incident  angle  can  be  any  desired  value  since  L/X  is  a  real 

number  and  not  limited  to  an  integer.  Additionally,  average  scattered  field  intensity 
is  independent  of  surface  length  as  long  as  many  correlation  lengths  are  included. 
Truncating  the  Floquet  expansion  in  equation  (2.2)  for  the  scattered  field  in  region  I 
and  expressing  it  in  Fourier  series  such  that  all  propagating  modes  and  2Ne  number 
of  evanescent  waves  are  included  results  in 


Ni 


Bm  -  m  $ ( 


(2.3) 


k’  = 


\A *-(¥)*. 


with 
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where  jV3  =  int(L/ A)  +  Ne  and  int(L/ A)  truncates  L/A  to  an  integer.  Writing 
equation  (2.3)  in  the  sine  and  cosine  form  of  a  Fourier  series  leads  to 

ill  i  f2jrm0  I  ,  .  .  f2ffm0  1 

ip  (x,z)  =  cos  I  x  -  +  ism  I — — — x  -  k;.z\  +  a0e  '  ’ 

+  2  [a™ cos  (~p*)  +  bm  sin  (“T"1)]  e'k‘m(‘~i>  (2-4) 

For  region  II  (d  >  z  >  f(x)),  the  finite  element  method  is  applied  to  determine  the 
modal  fields  4ij^(x,z)  where  m  is  the  modal  index.  Modal  fields  obey  the  following 
equations: 


(v2  +  £2)tfc(x,a)  =  0,  m  =  0,1,2,...,%  (2.5) 

(v*  +  Jfc2) #?••(*,*)  =  o,  mol, 2,3,...,%  (2.6) 

where  superscripts  c  and  s  stand  for  cosine  and  sine  respectively.  Respective  boundary 


conditions  for  each  are: 

dWixJjx))  Q  (27) 

e(-i/V)  =  tfW2,»)  (2-8) 

<0  =  cos  ("^"*)  (2.9) 

wt&m=0  (2.io) 

^(-L/2,z)  =  ^(L/2,z)  (2.11) 

<l’m’(x,d)  =  sm(^—x'j  (2.12) 


Equations  (2.7)  and  (2.10)  specify  the  Neumann  boundary  conditions  pnd  distinguish 
the  work  accomplished  in  this  thesis  from  that  of  references  [8]  and  (9). 

The  wave  equations  (2.5)  and  (2.6)  along  with  equations  (2.7)-(2.12)  allow  conve¬ 
nient  solution  for  and  by  the  finite  clement  method  for  each  modal  index 
m.  Discretization  of  region  II  is  finite,  and  the  modal  field  solutions  are  real,  thus, 
only  a  real,  sparse  matrix  solver  is  required  (18). 
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After  solving  for  the  unknown  modal  fields,  the  total  field  in  region  II  can  be 
expressed  as: 

Wj 

*"(*,*)  =  F0^(x,z)  +  £  [Fmtf*(x,z)  +  <?„tf •*(*,*)  (2.13) 

m=l 


allowing  the  real  boundary  conditions  of  equations  (2.9)  and  (2.12)  to  be  used. 

The  required  number  of  equations  to  determine  the  unknown  modal  amplitudes  Bm 
are  obtained  by  matching  the  continuities  of  ip1  and  tp11  and  their  normal  derivatives 
at  z  =  d.  From  equations  (2.4)  and  (2.13)  and  the  orthogonal  properties  of  the  sine 
and  cosine  functions,  the  following  equations  are  obtained: 


A  e-,b 

f,<i  +  Or 

n  —  — 0,1,.. 

;N3 

(2.14) 

M  e~'b 

^  +  5„ 

i  ~  Gjny  771  =  1,2,.. 

■  ,n3 

(2.15) 

i/1 

LJa 

dxFo^-—-(x,d) 

(2.16) 

1  ti 

VtZL 

"  M-l  *'0 


dx  \Fm 


dz 


(x,d)  +  Ga;-~ - 


(x,d) 


,  OPo* 


+  =  j  J  dxcos  (-J-*)  F0^—-(x,d) 


(2.17) 


/:-“(¥■) 

x[FM^p(x,d)  +  GM^(x,d) 


,  m  0 


'*  4*  F'.m — 


.  /'2xm  \ 

""IT 


2  /t-  ,  .  /2xm 
—  /  axsm  ( 

£  JO 


F0 


(x,d) 


(2.18) 


2  £2,  /t,  _ 

+z  sy#  * 

,ii* 


&) 


r.  _  dlpM*  . 

Fn-g^-(x,d)  +  Gjii—g^-(x,d) 


,  m  ^  0 
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where  Smmc  is  the  Kronecker  delta.  Substituting  equations  (2.14)  and  (2.15)  into 
(2.16),  (2.17),  and  (2.18)  leads  to  the  following  equations  in  terms  of  the  unknowns 
^0)  fim;  and  bm. 

— *2 kao  +  a0—  J  dx^~—{x,  d) 

+ 1,  {*“  [U^-- f^]  +  i"  [iJo1’- f^)]} 

e  ,k,,i l^-i2k,t6omo  —  Som<,-£ Jo  dx^~-(x,d)  (2.19) 

-ik'mam  +  fl0|j[  dx cos  ~p(x,d) 

+sMy^dxcoi^x)d-f-{x’d\ 

~  e  '*"*  j-iMmm, dxcos^^xj  ^2— (X)<f)  (2.20) 

2  /l  ,  f2xm  \  d4>!J‘c 

~lI  dxcm{-r*)-it‘MU*«-i 
~ilJodxCOS  (" rx )  ,  m  ^  0 

->Om  +  «of^fasin  (^x)  ^(x,d) 

=  -  SQmlj^dxsm  (^V)  ^(x,d)  (2.21) 
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where  Umo-i  —  0  for  normal  incidence  (mo  =  0)  and  Um j_i  =  1  for  oblique  incidence 
(mo  5^  0).  In  equations  (2.19),  (2.20),  and  (2.21),  ^g-(x,d)  and  —£‘-(x,d)  are 
evaluated  using  an  IMSL  subroutine  which  takes  the  derivatives  of  quadratic  inter¬ 
polation  of  $[*’'(*,  z)  and  z)  along  the  z  direction.  Integrations  are  evaluated 

using  Simpson’s  rule.  After  a o,  am,  and  bm  are  determined,  the  unknown  modal  field 
amplitudes,  Bm,  can  be  determined  by  the  following: 


B-mt}  —  no 


(2.22) 


Bm.m  =  l(am  -  ibm),  N3>m>0  (2.23) 

B-m- mo  =  ^(am  +  ibm),  N3  >  m  >  0  (2.24) 


2.2  Implementation 

Neumann  boundary  conditions  are  a  natural  result  of  the  Galerkin  procedure  applied 
to  the  finite  element  method  here.  In  the  finite  element  implementation,  three  general 
factors  affect  the  accuracy  of  the  solution:  (i)  Nz,  the  number  of  nodes  along  the 
x-axis,  (ii)  A’,,  the  number  of  nodes  along  the  z-axis,  and  (iii)  2Nb,  the  number  of 
evanescent  modes  considered.  The  parameter  d  must  be  set  to  be  above  the  maximum 
height  of  the  rough  surface  allowing  some  variation  of  the  field  from  its  surface  value 
to  its  value  at  z  =  d,  yet  small  enough  to  minimize  the  number  nodes  in  the  z 
direction.  Convergence  of  the  method  is  verified  by  increasing  Arx,  A,,  and  2Nb-  As 
it  turns  out,  the  number  of  evanescent  modes  is  not  a  sensitive  parameter  relative 
to  power-conservation  tests  (9)  For  the  numerical  examples  examined  with  the  finite 
element  method,  the  following  parameter  values  are  used,  2Ne  =  4,  N,/X  w  10  and 
Nx/I  zt  7  for  l  =  0.4A  and  Nz/X  zs  10  for  /  =  A.  The  parameter  d  is  set  at  0.1A  above 
the  maximum  surface  height  for  each  realization. 

Following  standard  finite  clement  procedure,  region  II  is  divided  into  a  number 
of  first-order,  triangular  elements.  The  modal  field  for  each  mode  m  within  each 
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element  (e)  is  expressed  in  terms  of  the  field  at  the  three  nodes  of  each  triangle  in 
the  following  matrix  notation: 

1&W)  =  [A'WJ  {^M}  (2.25) 

Field  values  at  the  three  nodes  are  represented  by  |jVMj  are  linear  interpolation 
functions  for  each  element  (e),  and  {•}  and  [-J  denote  a  column  and  row  vector 
respectively.  Complete  representation  of  the  modal  field  in  region  II  consisting  of  M 
elements  follows: 

M  M 

*)  =  £  &«W)  =  £IA'WJ  {4>(c)}  (2-26) 

e=l  e=l 

After  substituting  equation  (2.26)  into  equations  (2.5)  and  (2.6)  and  applying  the 
Galerkin  procedure,  the  following  results: 

MW  =  {/»}  (2.27) 


where 


Ml 


c WW  1  fcW<c>)  I  dN^  I 
dx  J  \  dz  )  [  dz  J 

-k2  J  ^{iV'}[A'"Jdidcj 


dxdz 


Region  II  has  a  total  of  NzNt  free  and  prescribed  nodes  and  2(Arr  -  1)(AI»  -  1) 
number  of  elements.  There  are  Nx  number  of  prescribed  nodes  at  the  boundary  z  =  d. 
Neumann  boundary  conditions  are  a  natural  result  of  the  Galerkin  procedure.  Hence, 
surface  nodes  are  considered  free  nodes  and  surface  node  field  values  are  determined 
in  the  same  manner  as  nodal  field  values  in  the  interior  of  region  II.  For  periodic 
boundary  conditions  at  x  and  x  +  L,  field  values  must  be  identical  but  their  values 
are  unknown.  The  periodicity  is  satisfied  by  setting  the  global  node  numbers  at  x 
and  i+i  to  the  same  value  for  a  fixed  z.  As  a  result,  ^  in  (2.29)  is  a  (A'x— l)(Afx  — 1) 
column  vector  of  free  nodes  only,  [A]  is  a  (Arx-1)(A'I-1)  by  (AX-1)(A’X-1)  real  and 
symmetric  matrix  with  half-bandwidth  2Nt  corresponding  to  the  global  assembly  of 
free  nodes.  P  is  a  (AC  —  l)(Afx  —  1)  column  vector  of  assembled  prescribed  nodes. 
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Equation  (2.27)  must  be  solved  once  for  each  of  the  (2 iV3  +  1)  modes  in  (2.5)  and 
(2.6).  A  direct  sparse  matrix  solver  utilizing  LU  decomposition  and  backsubstitution 
is  used  (18).  Since  [A]  is  sparse  and  symmetric,  profile  storage  of  its  upper  half  only 
significantly  reduces  memory  storage  requirements.  The  square  matrix  [A]  contains 
only  free  nodes  and  remains  unchanged  for  all  of  the  modes,  hence,  the  matrix  solver 
accomplishes  LU  decomposition  once  and  backsubstitution  (2Ar3  +  1)  times. 

Results  of  the  TWIE  method  are  computed  in  terms  of  the  bistatic  scattering 
cross  section  as  defined  in  (16).  For  convenience,  TWIE  results  are  plotted  using  a 
normalized  cross  section  o(6,)/  cos  0,  so  that  the  integral  over  all  scattered  angles  0, 
reduces  to  unity  for  power  conservation.  In  the  discrete  case  of  the  finite  clement 
method,  the  sum  of  scattered  power  over  all  discrete  scattered  angles  0m  must  re¬ 
duce  to  unity.  Since  2m\Bm\2  cosOm/  cosB,  =  / dBmLcosOm\Bm\^  cos  0m/cos0i,  the 
normalized  bistatic  scattering  cross  section  for  all  discrete  scattered  angles  is  defined 
as  L\Bm  |2  cos2  0m/ cos  0i. 

2.3  Decomposition  into  Coherent  and  Incoherent  Scattering 

Analysis  of  the  effects  of  decreasing  surface  length  on  finite  element  results  over  100 
realizations  is  accomplished  by  plotting  the  incoherent  scattering  cross  section.  This 
is  an  important  result  since  it  is  the  measured  quantity  in  monostatic  and  bistatic 
scattering.  The  incoherent  scattering  cross  section  is  calculated  by  first  determining 
the  incoherent  intensity  as  defined  in  the  following: 

<  l\&m|2  >  -I  <  i>m  >  ?  (2.28) 

where 

=  Bmeik‘",*lk‘"^  (2.29) 

<  ItM2  >  is  the  total  scattered  intensity,  |  <  ipm  >,|2  is  the  coherent  scattered 
intensity,  and  m  =  0, ±1,±2, ....  Since  the  parameter  d  varies  with  each  surface 
realization  depending  upon  its  respective  maximum  height,  the  phase  term 

must  be  considered  in  the  incoherent  scattering  cross  section  calculation.  Therefore, 
after  the  modal  scattered  fields  Bm  are  determined  for  each  surface  realization,  their 
corresponding  tjrm  values  are  determined  at  z  =  dm„.  The  value  for  dmaz  is  defined 
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as  the  maximum  d  value  for  all  100  realizations.  Hence,  each  field  value  and  the 
subsequent  incoherent  scattering  cross  section  calculation  is  referenced  to  a  common 
phase  point  at  the  maximum  d. 


Chapter  3 

PERTURBATION  THEORY 

The  Rayleigh-Rice  small  perturbation  theory  (SP)  is  typically  thought  to  be  valid 
for  slightly  rough  surfaces.  More  specifically,  the  rms  surface  height  h  should  be  small 
compared  to  wavelength  of  the  incident  field.  This  is  usually  expressed  as  kh  <.  1 
with  k  being  the  wavenumber  [17].  A  recent  study  by  Thorsos  in  reference  [17]  better 
clarifies  the  limits  of  validity  of  the  first  order  perturbation  theory  which  is  used  here 
for  comparison  with  the  finite  element  method. 

3.1  Formulation 

In  the  classical  approach  to  the  perturbation  theory,  the  scattered  field  is  represented 
as  a  superposition  of  outgoing  plane  waves  [17].  The  scattered  intensity  is  averaged  us¬ 
ing  the  statistical  properties  of  the  rough  surfaces.  Considering  plane  wave  incidence 
and  expressing  the  scattered  field  as  a  superposition  of  plane  waves,  the  following 
expression  for  the  scattered  field  in  the  region  above  the  rough  surface  results: 

*,=  £  Bmeik”*+ik"‘  (3.1) 

m=-oo 

The  unknown  scattered  plane  wave  amplitudes  are  represented  by  Bm.  Applying 
Neumann  boundary  conditions  on  the  surface,  expanding  in  smallness,  and  averag¬ 
ing  according  to  statistical  theory,  the  following  expression  results  for  the  average 
scattered  intensity: 
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L  is  the  surface  length,  h  is  the  rms  surface  height,  1  is  the  correlation  length,  and 
mo  is  the  discrete  incident  angle  index  as  in  the  finite  element  analysis.  Both  the 
coherent  and  incoherent  contributions  to  the  average  scattered  intensity  are  included 
in  equation  (3.2)  with  the  Kronecker  delta  function  representing  the  coherent  com¬ 
ponent.  The  one-dimensional  equivalent  to  the  normalized  bistatic  scattering  cross 
section  required  for  comparison  with  finite  element  results  follows  from  the  average 
scattered  intensity: 


o{0m) 


Lcos20m  f ,  .  .[hi  -  kaxk„  -  i] 

-^r\Smm°+A — tj - —e  j 


(3.3) 


Perturbation  results  are  presented  in  terms  of  discrete  scattered  angles  in  an  identical 
manner  to  the  finite  element  method. 


I 
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Chapter  4 

KIRCHOFF  APPROXIMATION 


The  Kirchoff  approximation  (KA)  is  generally  thought  to  apply  to  gently  undu¬ 
lating  sufaces  which  are  defined  as  surfaces  with  a  radius  of  curvature  that  is  large 
compared  to  a  wavelength.  This  criterion  places  no  general  restrictions  upon  the 
surface  height,  and  the  Kirchoff  approximation  has  been  applied  to  higher  frequency 
cases  than  the  perturbation  theory  as  a  result.  Recent  investigations  by  Thorsos 
better  quantify  the  domain  of  validity  for  the  Kirchoff  approximation  as  presented  in 
reference  (16]. 


4.1  Formulation 


The  classical  starting  point  for  the  Kirchoff  approximation  is  the  Helmholtz  integral 
which  defines  the  scattered  field  in  relation  to  the  surface  field  as  presented  in  the 
following: 

_,'dGo(r,r')  ,  r')] 


VH1-)  =  J ds  ^(r')dgjv—  -  G0( r,  r')- 


dn‘ 


('1-1) 


The  following  substitutions  are  made  for  the  appropriate  far  field  Green’s  function, 
Neumann  boundary  conditions,  and  assumed  plane  wave  incidence: 


Go  i 


(4.2) 


n 

dn 


=  e’kir' 


(4.3) 

(4.4) 

(4.5) 


Theoretical  averaging  of  the  scattered  intensity  is  accomplished  using  the  statistical 
properties  of  the  rough  surface  as  in  the  perturbation  theory.  Combining  the  coherent 
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and  incoherent  contributions  to  the  normalized  bistatic  scattering  cross  section,  the 
following  results: 


°(Pm) 


n=l 


M)2n  -Li 

n\y/n  &  " 


where 


,  2  irm0 

“  -  L 

v  =  &(sind;  —  sin  $m)x— &(««£>.•  + cos  0m)z 
ks  =  ksmOmii  +  fccos0mz 
vx  =  k( sin  0,  —  sin  6m) 
vz  —  k(cos  6i  +  cos  6m) 


(4.6) 


Results  are  presented  in  terms  of  discrete  scattered  angles  in  an  identical  manner  to 
the  finite  element  method. 


Chapter  5 

NUMERICAL  RESULTS 


Numerical  results  and  comparisons  of  FEM,  TWIE,  SP,  and  KA  are  presented 
in  this  chapter.  Section  5.1  provides  convergence  test  results  for  FEM  with  respect 
to  the  number  of  nodes  in  both  the  x  and  z  directions.  Section  5.2  compares  the 
FEM  and  TWIE  result  for  the  analysis  of  a  single  surface  realization.  Comparisons 
between  FEM  and  TWIE  with  respect  to  required  computer  CPU  time  and  memory 
storage  are  presented  in  section  5.3.  Section  5.4  presents  a  comparison  of  each  of  the 
four  approaches  for  various  surface  characteristics.  Power  conservation  test  results 
are  compared  for  FEM  and  TWIE  in  section  5.5.  Finally,  section  5.6  presents  an 
analysis  of  the  effects  of  decreasing  surface  length  on  FEM  results. 

For  each  FEM  and  TWIE  result  in  sections  5.4,  5.5,  and  5.6,  100  surface  realiza¬ 
tions  are  generated  with  predetermined  rms  surface  height,  A/A,  correlation  length, 
I/A,  and  surface  length,  A/A.  The  rms  surface  slope  s,  rms  height,  correlation  length, 
and  rms  slope  angle  7  are  related  by  s  =  \/2h/l  =  tan  7.  The  FEM  and  TWIE 
programs  are  then  executed  in  order  to  calculate  the  ensemble  average  over  the  real¬ 
izations.  In  the  TWIE  method,  the  incident  wave  consists  of  a  plane  wave  which  is 
tapered  using  a  Gaussian  taper  function.  The  result  is  an  incident  wave  consisting 
of  an  angular  spectrum  of  plane  waves  about  the  incident  angle  with  a  width  of  ap¬ 
proximately  4A/(v/2jrAcos0,)  (16).  TWIE  results  presented  in  this  thesis  are  based 
on  the  computer  code  of  reference  (4). 

5.1  Convergence  of  FEM  with  Respect  to  Nx  and  Nz 

Figure  5.1  illustrates  the  effect  on  the  finite  element  result  of  increasing  the  number 
of  nodes  in  the  x  direction  from  Nx  =  301  to  Nz  =  601  for  a  single  realization  with 
0;  =  0°,  A/A  =  0.1,  I/A  -  1.0,  and  A/A-  =  30.2.  Figure  5.2  illustrates  the  effect  of 
increasing  the  number  of  nodes  in  the  z  direction  from  N,  =  8  to  IV,  =  15  for  the 
same  realization.  In  each  case,  there  is  virtually  no  noticeable  difference  between  the 
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respective  dB  plots.  Hence,  these  two  plots  demonstrate  convergence  of  the  FEM 
result  with  respect  to  the  number  of  nodes  in  both  the  x  and  z  directions. 


5.2  Comparison  of  FEM  and  TWIE  for  a  Single  Realization 

Figure  5.3  compares  the  result  of  FEM  and  TWIE  analysis  of  the  same  surface  re¬ 
alization  used  to  test  FEM  convergence  with  respect  to  increasing  Nz  and  A*.  It  is 
evident  that  differences  do  exist  between  the  two  methods  for  a  single  realization. 
This  is  to  be  expected  as  the  FEM  and  TWIE  approaches  effectively  sec  two  different 
realizations  since  the  former  has  periodic  boundary  conditions  and  the  latter  utilizes 
a  tapered  incident  wave.  With  surface  length  L/ A  =  30.2,  the  incident  wave  has  an 
angular  width  of  about  ±1.5°  with  0,  =  0°  for  TWIE.  Despite  differences  between 
FEM  and  TWIE  for  a  single  realization,  results  are  excellent  upon  averaging  over 
many  realizations  as  illustrated  in  section  5.4. 


5.3  Comparisons  of  CPU  Time  and  Memory  Storage 

In  Figures  5.4  and  5.5,  CPU  time  and  memory  storage  requirements  for  FEM  and 
TWIE  approaches  are  compared.  Results  are  based  upon  the  use  of  a  VAX  Station 
3500  for  a  single  realization  with  h/X  =  0.1  and  l/X  =  0.4.  Plots  are  made  as  a 
function  rf  surface  length  with  CPU  time  presented  in  seconds  and  memory  storage 
presented  in  Megabytes.  It  is  evident  from  the  plots  that  TWIE  generally  requires 
more  CPU  time  and  memory  storage  than  FEM.  This  is  true  when  considering  scat¬ 
tering  problems  where  dimensions  of  the  scattering  geometry  are  much  greater  in  one 
direction  than  in  the  other.  Scattering  from  a  one-dimensional  surface  with  large 
surface  length  meets  this  criterion,  hence,  the  sparsity  of  the  real  FEM  matrix  for 
solving  equation  (2.27)  results  in  significant  savings  in  CPU  time  and  memory  stor¬ 
age  as  compared  to  TWIE  which  requires  full,  complex  matrix  solution  and  storage. 
The  disparity  between  the  two  methods  becomes  more  significant  as  surface  length 
increases,  and  it  should  be  noted  that  a  large  surface  length  is  required  in  the  TWIE 
approach  in  order  to  minimize  the  angular  spectrum  of  the  incident  field  and  more 
accurately  predict  the  scattering  cross  section. 
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5.4  Comparison  over  100  Realizations 

Figure  5.6  compares  the  analytical  perturbation  theory  to  the  FEM  and  TWIE  results 
averaged  over  100  surface  realizations  with  5,-  —  44.06°  (m0  =  21),  h/X  =  0.05, 
l/X  =  0.4,  and  L/ X  =  30.2.  Surface  characteristics  are  in  the  domain  of  validity  of 
the  small  perturbation  theory,  and  results  agree  well  with  FEM  and  TWIE  results. 
A  comparison  between  the  small  perturbation  theory  and  FEM  for  Neumann  and 
Dirichlet  boundary  conditions  is  presented  in  Figure  5.7.  Dirichlet  results  are  based 
on  the  finite  element  code  of  reference  (9).  Surface  characteristics  are  the  same  as 
in  the  previous  illustration.  Neumann  boundary  conditions  result  in  a  more  diffuse 
scattered  field  than  that  of  the  Dirichlet  as  is  evident  in  Figure  5.7.  Results  in  Figure 
5.8  compare  the  Kirchoff  approximation  to  FEM  and  TWIE.  Surface  realizations  are 
such  that  0j  =  0°  (mo  =  0),  h/X  =  0.1,  l/X  =  1.0,  and  L/X  —  70.5.  The  Kirchoff 
approximation  is  considered  valid  in  this  case,  and  results  are  excellent  except  for 
scattered  grazing  angles  beyond  ±50°  where  I<A  underpredicts  the  scattered  intensity. 
Figure  5.9  illustrates  the  effect  of  increasing  the  rms  slope  angle  7  which  becomes 
19.47°  as  compared  to  8.05°  in  Figure  5.8.  As  expected,  the  angular  distribution  of  the 
incoherent  component  to  scattering  cross  section  is  increased.  With  0i  =  44.06°  (mo  = 
21),  h/X  =  0.1,  l/X  =  0.4,  and  L/X  =  30.2,  KA  results  are  excellent  for  0,  >  20°,  but 
I<A  significantly  underpredicts  scattered  intensity  for  0,  <  20°.  In  Figure  5.10,  h/X 
has  been  increased  from  0.1  to  0.2.  As  expected,  the  coherent  component  (specular 
peak)  is  significantly  reduced  from  that  of  Figure  5.9.  The  KA  results  are  generally 
poor  except  for  10°  <  0,  <  47°  although  improvement  over  Figure  5.9  is  evident 
for  0,  <  0°.  The  rms  height  h/X  is  increased  to  0.4  in  Figure  5.11  which  results  in 
7  =  54.74°.  KA  is  no  longer  valid,  hence,  only  FEM  and  TWIE  results  are  presented. 
As  in  all  previous  plots  averaged  over  100  realizations,  agreement  between  TWIE  and 
FEM  is  excellent.  The  specular  peak  has  completely  subsided  in  Figure  5.11  which 
is  not  surprising  given  the  rather  large  rms  height  specified. 

5.5  Comparison  of  Power  Conservation  Tests 

Table  5.1  compares  power  conservation  results  of  FEM  and  TWIE  for  Figure  5.6  and 
Figures  5.8-5.10.  Negative  values  indicate  an  ovcrprediction  of  scattered  power,  while 
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positive  values  indicate  an  underprediction.  In  general,  FEM  has  less  error  in  terms 
of  power  conservation  than  TWIE  with  the  largest  value  being  only  -7.9  x  10'3 
percent. 

5.6  FEM  Result  Dependence  on  Surface  Length 

In  order  to  investigate  the  effects  of  decreasing  surface  length  on  FEM  results,  further 
tests  with  100  surface  realizations  are  conducted.  The  coherent  component  to  the 
scattering  cross  section  is  subtracted  to  avoid  plotting  a  specular  peak  which  broadens 
with  decreasing  surface  length.  Incoherent  intensity,  <  |0m|2  >  — |  <  i>m  >  |2,  is 
calculated  at  the  maximum  d  for  all  100  surface  realizations.  Surface  realization 
characteristics  A/A  =  0.141  and  /  =  0.4  for  all  surface  lengths  are  chosen  to  produce 
a  relatively  high  rms  slope  of  s  =  0.5.  High  rms  slope  for  the  realizations  insures 
a  broad  angular  distribution  of  the  incoherent  result.  Figure  5.12  illustrates  the 
removal  of  the  coherent  contribution  to  the  total  scattering  cross  section  for  L/X  = 
30.2.  A  depression  in  the  incoherent  scattering  cross  section  at  0m  =  0;  results  after 
removing  the  specular  peak.  This  occurrence  is  unexplained  by  scattering  theory,  and 
more  work  is  required  to  determine  why  it  exists.  The  incoherent  result  of  Figure 
5.12  is  plotted  against  results  of  all  other  surface  lengths  tested.  Surface  length  is 
decreased  in  approximately  5A  intervals  from  L/X  =  30.2.  Actual  surface  length 
values  are  arrived  at  through  consideration  of  integer  values  for  mo  in  order  that 
0{  ss  44.06°  for  all  surfaces  tested.  Table  5.2  displays  the  specific  criteria  chosen 
for  each  test.  Figures  5.13-5.17  compare  the  incoherent  result  for  respective  surface 
lengths  of  24.45A,  20.13A,  14.38A,  10.07A,  and  4.31A.  An  IMSL  subroutine  is  used 
to  interpolate  points  for  L/X  <  14.38.  Excellent  agreement  can  be  seen  between 
results  down  to  L/X  =  10.07.  Results  for  L/X  =  4.31  are  also  excellent  for  0m  <  30°. 
However,  due  to  the  limited  number  of  scattered  angles  for  which  scattering  cross 
section  values  are  calculated,  the  dip  at  0m  =  Oi  has  a  significant  effect.  Table  5.2  also 
compares  required  computer  storage  and  CPU  time  for  each  test.  Memory  storage  is 
presented  in  megabytes  for  a  single  realization,  and  CPU  time  is  presented  in  seconds 
for  the  full  100  realizations  on  a  VAX  6000-440  workstation.  Significant  savings 
between  L/X  =  30.2  and  L/X  =  10.07  in  storage  and  CPU  time  is  evident.  Test  results 
are  encouraging  for  eventual  FEM  analysis  of  two-dimensional  surface  scattering. 
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Since  the  surface  length  required  for  accurate  FEM  results  is  relatively  small,  it  is 
possible  that  CPU  time  and  memory  storage  could  be  kept  within  reasonable  limits. 
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Figure  5.7:  Comparison  of  SP  and  FEM  for  Neumann  and  Dirichlet  Boundary 
Conditions:  (<?,•  =  44.06°,  h/X  =  0.05,  l/X  =  0.40,  LjX  =  30.20) 
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Table  5.2:  FEM  Comparison  for  Different  Surface  Lengths 


Figure  # 

L/X 

m0 

6, 

E9 

ea 

Memory  (MB) 

CPU  Time  (sec) 

5.12 

30.20 

m 

44.06° 

ex 

1.41 

1.108  x  104 

5.13 

24.45 

Q 

44.05° 

Esa 

El 

1.09 

6.750  x  103 

5.14 

20.13 

KOI 

44.07° 

m 

El 

0.87 

4.482  x  103 

5.15 

14.38 

m 

44.06° 

f3i 

EX 

0.59 

2.326  x  103 

5.16 

10.07 

D 

44.04° 

EX 

0.40 

1.285  x  103 

5.17 

4.31 

D 

44.11° 

m 

EX 

0.16 

3.561  x  103 
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Figure  5.12:  Incoherent  and  Total  Scattering  Cross  Sections:  (0; 
h/X  =  0.141, 1/X  =  0.40,  h/X  =  30.20) 
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Figure  5.14:  Comparison  of  FEM  for  L/\  =  30.20  and  L/X  =  20.13: 
//A  =  0.40) 
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Figure  5.16:  Comparison  of  FEM  for  L/X  =  30.20  and  Lj X  =  10.07:  (/i/A  —  0.141, 
1/ X  =  0.40) 


Chapter  6 

CONCLUSIONS 


The  finite  element  method  is  applied  to  rough  surface  scattering  with  periodic, 
Neumann  boundary  conditions.  Convergence  of  the  result  is  verified  by  varying  the 
number  of  nodes  in  both  the  x  and  z  directions.  Accuracy  is  proven  through  compar¬ 
ison  with  “exact”  results  of  an  integral  equation  solution  and  by  power  conservation 
tests.  Comparisons  with  the  Rayleigh-Rice  small  perturbation  theory  and  Kirchoff 
approximations  are  also  made  in  their  respective  regions  of  validity.  Finite  element 
advantages  over  the  integral  equation  approach  in  both  computer  CPU  time  and 
required  memory  storage  are  demonstrated. 

Additional  tests  show  that  the  surface  length  necessary  for  accurate  incoherent 
scattering  cross  section  for  the  finite  element  solution  is  as  little  as  10.07A.  As  a 
result,  substantial  savings  in  required  memory  storage  and  CPU  time  are  realized. 
Additional  research  is  needed  in  order  to  explain  the  resulting  depression  in  the  in¬ 
coherent  scattering  cross  section  at  0m  =  0;.  The  depression  which  existed  after 
removal  of  the  specular  peak  is  unexpected  and  is  not  explained  by  scattering  theory. 
However,  results  obtained  are  significant  as  they  make  it  likely  that  the  finite  ele¬ 
ment  method  will  prove  to  be  a  feasible  approach  to  large-scale  scattering  problems 
including  two-dimensional  surface  scattering. 
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